# -*- coding: utf-8 -*-
"""
Created on Thu May  5 15:08:00 2022

@author: 17136
"""

#无缝线路检算
import xlrd
import numpy
import math
data=[]
#读取参数表格
wb=xlrd.open_workbook(r'C:\Users\10707\PycharmProjects\pythonProject\2-canshu220231003-q.xlsx')
sheet=wb.sheets()[0]
for a in range(1,sheet.nrows):
    cells=sheet.row_values(a)
## 输入参数
    CZ=cells[0]#车站名称
    GH=cells[1]   #股道编号
    JC=cells[2]#机车类型
    Gd=cells[3]   #轨道类型
    GZ=cells[4]#轨枕类型
    ZJ=cells[5]   #轨枕级数
    GG=cells[6]#钢轨类型
    GM=cells[7]   #钢轨磨耗
    GL=cells[8]#钢轨种类
    Kn=cells[9]   #轨下胶垫刚度
    Tx=cells[10]#锁定轨温修正值（取0）
    h=cells[11]   #既有道床厚度
    R=cells[12]#曲线半径
    v=cells[13]   #线路允许速度
    ch=cells[14]#设计超高
    Tmax=cells[15]   #当地最高轨温
    Tmin=cells[16]#当地最低轨温
    Hzf=cells[17]   #荷载方式
    P1=cells[18]#一位轮重
    P2=cells[19]   #二位轮重
    P3=cells[20]#三位轮重
    Pa1=cells[21]   #轮位间距1
    Pa2=cells[22]#轮位间距2
    Te=cells[23]   #设计锁定轨温
#

## 计算参数确定
    if GZ==2:
        L=2500;     #轨枕长度
        a1=500;     #荷载作用点至轨端距离
        b=262.5;      #轨枕底面宽度
        Hzl=7.0*ZJ/1760;    #线路纵向阻力值
        if Kn==110:         #轨下胶垫刚度
            D=30000;        #钢轨支撑刚度
        else:
            D=27200;

        if R<=800:
            Q=8.9*(ZJ/1840)+8.9*(ZJ/1840)*0.1;  ###
        else:
            Q=8.9*(ZJ/1840);

    if GZ==22:
        L=2500;     #轨枕长度
        a1=500;     #荷载作用点至轨端距离
        b=265;      #轨枕底面宽度
        Hzl=8.8*ZJ/1760;    #线路纵向阻力值
        if Kn==110:         #轨下胶垫刚度
            D=30000;        #钢轨支撑刚度
        else:
            D=27200;

        if R<=800:
            Q=8.9*(ZJ/1840)+8.9*(ZJ/1840)*0.1;
        else:
            Q=8.9*(ZJ/1840);

    if GZ==3:
        L=2600
        a1=550
        b=290
        D=33000;
        Hzl=15*ZJ/1667;
        if R<=800:
            Q=11.5*(ZJ/1667);
        else:
            Q=11.5*(ZJ/1667);

    a=1000000/ZJ;
    Pj=490000;
    e=950;
    e1=3*L/8+e/4;
    M=1.6;#道床不均匀系数
    gal=0.15;#路基面允许应力
    gad=0.5;#道床允许应力
    Dk=40;#道床肩宽
    Hk=150;#砟肩堆高
    if GG==50:
        dteP=130;
        b1=132;
        ag=18;
        if GM==0:
            F=6580;
            W1=250000;
            W2=289000;
            Ix=20370000;
            Iy=3770000;
        elif GM==3:
            F=6580;
            W1=242000;
            W2=283000;
            Ix=19460000;
            Iy=3770000;
        elif GM==6:
            F=6580;
            W1=230000;
            W2=275000;
            Ix=18270000;
            Iy=3770000;
        elif GM==9:
            F=6580;
            W1=216000;
            W2=264000;
            Ix=17020000;
            Iy=3770000;
    elif GG==60:
        detP=154;
        b1=150;
        ag=18;
        if GM==0:
            F=7745;
            W1=339400;
            W2=396000;
            Ix=32170000;
            Iy=5240000;
        elif GM==3:
            F=7745;
            W1=318000;
            W2=385000;
            Ix=30690000;
            Iy=5240000;
        elif GM==6:
            F=7745;
            W1=291000;
            W2=375000;
            Ix=28790000;
            Iy=5240000;
        elif GM==9:
            F=7745;
            W1=264000;
            W2=363000;
            Ix=26900000;
            Iy=5240000;
    elif GG==75:
        detP=188;
        b1=150;
        ag=20;
        if GM==0:
            F=9504;
            W1=432000;
            W2=509000;
            Ix=44890000;
            Iy=6610000;
        elif GM==3:
            F=9504;
            W1=420000;
            W2=496000;
            Ix=43280000;
            Iy=6610000;
        elif GM==6:
            F=9504;
            W1=405000;
            W2=482000;
            Ix=40890000;
            Iy=6610000;
        elif GM==9:
            F=9504;
            W1=390000;
            W2=480000;
            Ix=38980000;
            Iy=6610000;
    E=210000;
    arf=0.0000118;
    if GL==1:
        gag=457;
    else:
        gag=472;
    K=1.3;  #安全系数

    # 横向水平力系数
    if 300<=R<400:
        f0=-0.002*R+2.6;
    elif 400<=R<600:
        f0=-0.001*R+2.2;
    elif 600<=R<800:
        f0=-0.00075*R+2.05;
    elif 800<=R<2000:
        f0=1.45;
    elif 2000<=R<=20000:
        f0=1.3;
    elif R>20000:
        f0=1.25;
    #
    
    
    if JC==1:
        if v<=160:
            af0=0.6*v/100;
        else:
            af0=1;
    elif JC==2:
        if v<=160:
            af0=0.4*v/100;
        else:
            af0=1;
    bta=0.002*(11.8*v*v/R-ch);
    bta0=1;     #轨道框架刚度系数，有砟轨道可取1.0；
    l0=4000;
    R0=10/(1.39639*math.pow(10,-5));
    f=2;
    foe=3.575*l0*l0*math.pow(10,-8);
    fop=1.74549*l0*l0*math.pow(10,-7);
    t=foe/(l0*l0);
#轨道结构静力计算
    # 计算最大静位移、弯矩和枕上压力
    u1=D/a;
    u2=70000/a;
    gaz=10;
    k1=math.pow(u1/(4*E*Ix),1/4);
    k2=math.pow(u2/(4*E*Ix),1/4);
    y1=(k1/(2*u1))*1000*(P1*1+P2*math.exp(-k1*Pa1)*(math.cos(k1*Pa1)+math.sin(k1*Pa1))+P3*math.exp(-k1*(Pa1+Pa2))*(math.cos(k1*(Pa1+Pa2))+math.sin(k1*(Pa1+Pa2))));
    M1=(1/(4*k1))*1000*(P1*1+P2*math.exp(-k1*Pa1)*(math.cos(k1*Pa1)-math.sin(k1*Pa1))+P3*math.exp(-k1*(Pa1+Pa2))*(math.cos(k1*(Pa1+Pa2))-math.sin(k1*(Pa1+Pa2))));
    R1=((a*k2)/2)*1000*(P1*1+P2*math.exp(-k2*Pa1)*(math.cos(k2*Pa1)+math.sin(k2*Pa1))+P3*math.exp(-k2*(Pa1+Pa2))*(math.cos(k2*(Pa1+Pa2))+math.sin(k2*(Pa1+Pa2))));
    y2=(k1/(2*u1))*1000*(P1*math.exp(-k1*Pa1)*(math.cos(k1*Pa1)+math.sin(k1*Pa1))+P2*1+P3*math.exp(-k1*Pa2)*(math.cos(k1*Pa2)+math.sin(k1*Pa2)));
    M2=(1/(4*k1))*1000*(P1*math.exp(-k1*Pa1)*(math.cos(k1*Pa1)-math.sin(k1*Pa1))+P2*1+P3*math.exp(-k1*Pa2)*(math.cos(k1*Pa2)-math.sin(k1*Pa2)));
    R2=((a*k2)/2)*1000*(P1*math.exp(-k2*Pa1)*(math.cos(k2*Pa1)+math.sin(k2*Pa1))+P2*1+P3*math.exp(-k2*Pa2)*(math.cos(k2*Pa2)+math.sin(k2*Pa2)));             
    Mj=max(M1,M2);  #枕上弯矩
    Rj=max(R1,R2);  #枕上压力
    Md=Mj*(1+af0+bta)*f0;#钢轨最大动弯矩
    Rd=Rj*(1+af0+bta);  #轨枕最大动压力              
    #轨枕和道床厚度检算
    Mgd=Rd*(a1*a1/(2*e)-b1/8)*math.pow(10,-6);#轨下截面正弯矩Mg
    Med=-(Rd*(3*L*L+4*e*e-8*a1*e-12*a1*L)*math.pow(10,-6))/(4*(3*L+2*e));#枕中截面负弯矩
    gadc=Rd*M/(b*e1);#作用在道床的应力σ道床
    hdc=Rd/(2*e1*gal*math.tan(35*math.pi/180));#计算道床厚度h
    
 #无缝线路稳定性计算       
    #统一无缝线路稳定性计算公式
    Rx=R*1000*R0/(R*1000+R0);
    w=2*bta0*E*Iy*math.pi*math.pi*(t+4/(math.pi*math.pi*math.pi*Rx));
    ll=(w+math.sqrt(w*w+(4*Q/(math.pow(math.pi,3))-w*t/f)*2*f*bta0*E*Iy*math.pi*math.pi))/(4*Q/(math.pow(math.pi,3))-w*t/f);
    n=1;
    foe1=foe;
    while n<=10:
        foe1=foe1*ll/(l0*l0);
        t1=foe1/(l0*l0);
        w1=2*bta0*E*Iy*math.pi*math.pi*(t1+4/(math.pi*math.pi*math.pi*Rx));
        lb=(w1+math.sqrt(w1*w1+(4*Q/(math.pow(math.pi,3))-w1*t1/f)*2*f*bta0*E*Iy*math.pi*math.pi))/(4*Q/(math.pow(math.pi,3))-w1*t1/f);
        aa=ll;
        bb=lb;
        l0=math.sqrt(aa);
        ll=bb;
        n=n+1;
    Pw1=2*bta0*E*Iy*math.pi*math.pi*(f+foe1)/lb+4*Q*lb/(math.pow(math.pi, 3));
    Pw2=f+foe1+4*lb/(math.pi*math.pi*math.pi*Rx);
    Pw=Pw1/(Pw2*1000);
    Ppa=Pw/1.3;                     #最后得允许最大温度力
    Wys=Ppa*1000/(2*E*arf*F);       #轨道稳定性允许温升[ΔTu]
    
#无缝线路钢轨强度检算
    gat1=Md/W1;#轨头边缘最大可能动弯应力
    gat2=Md/W2; #轨低边缘最大可能动弯应力 
    Wyj=(gag/K-gat2-gaz)/(E*arf);   #轨道强度允许温降[ΔTd]
    Tu =Te + 3 + Tx;
    Td=Te - 3 - Tx;
    Tumax=Tmax-Td;
    Tdmax=Tu-Tmin;               
    gat=E*arf*Tdmax;  #钢轨温度拉应力σt
    gah=round(gat2+gat+gaz,2);#作用在钢轨的应力
    

#缓冲区长度计算：
    Ptmax=E*arf*Tdmax*F;#钢轨最大温度拉力
    Pcmax=E*arf*Tumax*F;#钢轨最大温度压力
    Lss=(Ptmax-Pj)/(Hzl*1000);#伸缩区长度Ls
    Lcsc=(Ptmax-Pj)*(Ptmax-Pj)/(2*E*F*Hzl);
    Lcsd=(Pcmax-Pj)*(Pcmax-Pj)/(2*E*F*Hzl);
    Lhsc=((Ptmax-Pj)*25000-0.25*Hzl*25000*25000)/(2*E*F);
    Lhsd=((Pcmax-Pj)*25000-0.25*Hzl*25000*25000)/(2*E*F);
    asx=ag-(Lcsd+Lhsd);#预留轨缝上限（冬季）a上
    axx=(Lcsc+Lhsc);#预留轨缝下限（夏季）a下
    ayl=(asx+axx)/2;#缓冲区预留轨缝长度a
    df=E*F*(arf*Tdmax)*(arf*Tdmax)/Hzl;#钢轨断缝λ
    data.append([CZ,GH,Te,gah,round(gag/K,2),round(Tumax,2),round(Wys,2),round(df,2),
                 round(Mgd,2),round(Med,2),round(gadc,3),round(hdc,1),h,round(Lss,1),
                 round(ayl,1),round(Ppa,1),round(Mj/1000000,2),round(Rj/1000,2)])

import xlwt
sc=xlwt.Workbook("jieguo")
sheet_jieguo=sc.add_sheet('sheet1',cell_overwrite_ok=True)  
title=['车站','股道','锁定轨温/℃','钢轨应力/MPa','评价结果','稳定性','评价结果','钢轨断缝/mm','评价结果','轨枕弯矩','评价结果','道床应力/MPa','评价结果','计算道床厚度/mm','评价结果','伸缩区长/m','预留轨缝/mm','允许温度力/kN','钢轨动弯矩/kN*m','轨枕动压力/kN']            
for k in range(len(title)):
    sheet_jieguo.write(0,k,title[k])
i=1
for m in data:
    sheet_jieguo.write(i,0,m[0])
    sheet_jieguo.write(i,1,m[1])
    sheet_jieguo.write(i,2,'%.1f±5' %m[2])
    if float(m[3])<float(m[4]):
        sheet_jieguo.write(i,3,'%.1f<%.1f' %(m[3],m[4]))
        sheet_jieguo.write(i,4,'满足')
    else:
        sheet_jieguo.write(i,3,'%.1f>%.1f' %(m[3],m[4]))
        sheet_jieguo.write(i,4,'不满足')
    if float(m[5])<float(m[6]):
        sheet_jieguo.write(i,5,'%.2f<%.2f' %(m[5],m[6]))
        sheet_jieguo.write(i,6,'满足')
    else:
        sheet_jieguo.write(i,5,'%.2f>%.2f' %(m[5],m[6]))
        sheet_jieguo.write(i,6,'不满足')
    if float(m[7])<=70:
        sheet_jieguo.write(i,7,'%.2f≤70' %m[7])
        sheet_jieguo.write(i,8,'满足')
    elif 70<float(m[7])<=90:
        sheet_jieguo.write(i,7,'%.2f>70' %m[7])
        sheet_jieguo.write(i,8,'超限')
    else:
        sheet_jieguo.write(i,7,'%.2f>90' %m[7])
        sheet_jieguo.write(i,8,'严重超限')
    if float(m[8])<=13.3 and abs(float(m[9]))<=10.5:
        sheet_jieguo.write(i,9,'%.1f≤13.3  |%.1f|≤10.5' %(m[8],m[9]))
        sheet_jieguo.write(i,10,'满足')
    elif float(m[8])>13.3 and abs(float(m[9]))<=10.5:
        sheet_jieguo.write(i,9,'%.1f>13.3  |%.1f|≤10.5' %(m[8],m[9]))
        sheet_jieguo.write(i,10,'不满足')
    elif float(m[8])<=13.3 and abs(float(m[9]))>10.5:
        sheet_jieguo.write(i,9,'%.1f≤13.3  |%.1f|>10.5' %(m[8],m[9]))
        sheet_jieguo.write(i,10,'不满足')
    else:
        sheet_jieguo.write(i,9,'%.1f>13.3  |%.1f|>10.5' %(m[8],m[9]))
        sheet_jieguo.write(i,10,'不满足')
    if float(m[10])<=0.5:
        sheet_jieguo.write(i,11,'%.2f≤0.5' %m[10])
        sheet_jieguo.write(i,12,'满足')
    else:
        sheet_jieguo.write(i,11,'%.2f>0.5' %m[10])
        sheet_jieguo.write(i,12,'不满足')
    if float(m[11])<=float(m[12]):
        sheet_jieguo.write(i,13,'%.1f≤%.1f' %(m[11],m[12]))
        sheet_jieguo.write(i,14,'满足')
    else:
        sheet_jieguo.write(i,13,'%.1f>%.1f' %(m[11],m[12]))
        sheet_jieguo.write(i,14,'不满足')
    if float(m[13])<=25:
        sheet_jieguo.write(i,15,25)
    elif 25<float(m[13])<=50:
        sheet_jieguo.write(i,15,50)
    elif 50<float(m[13])<=75:
        sheet_jieguo.write(i,15,75)
    elif 75<float(m[13])<=100:
        sheet_jieguo.write(i,15,100)
    sheet_jieguo.write(i,16,m[14])
    sheet_jieguo.write(i,17,m[15])
    sheet_jieguo.write(i,18,m[16])
    sheet_jieguo.write(i,19,m[17])
    i=i+1
sc.save('jieguo.xls')